Cartography with R - Exercises

Create reproducible maps with R
Published

March 11, 2025

Vector Data & Geoprocessing

Vector Data

sf package

Install the sf package.

install.packages("sf")

Import and export

Import municipalities.

library(sf)
#> Warning in CPL_gdal_init(): GDAL Message 1: GDAL AVIF driver was built against
#> libavif 1.1.1 but is running against 1.2.0. Runtime issues could occur
#> Linking to GEOS 3.13.1, GDAL 3.10.2, PROJ 9.5.1; sf_use_s2() is TRUE
#> WARNING: different compile-time and runtime versions for GEOS found:
#> Linked against: 3.13.1-CAPI-1.19.2 compiled against: 3.13.0-CAPI-1.19.0
#> It is probably a good idea to reinstall sf (and maybe lwgeom too)
st_layers(dsn = "data/lot.gpkg")
#> Driver: GPKG 
#> Available layers:
#>     layer_name geometry_type features fields              crs_name
#> 1     communes Multi Polygon      313     12 RGF93 v1 / Lambert-93
#> 2 departements Multi Polygon       96      5 RGF93 v1 / Lambert-93
#> 3  restaurants         Point      694      2 RGF93 v1 / Lambert-93
#> 4   elevations         Point     5228      1 RGF93 v1 / Lambert-93
#> 5       routes   Line String     1054      2 RGF93 v1 / Lambert-93
mun <- st_read("data/lot.gpkg", layer = "communes")
#> Reading layer `communes' from data source 
#>   `/home/tim/Documents/prj/cartography_with_r/data/lot.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 313 features and 12 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 539668.5 ymin: 6346290 xmax: 637380.9 ymax: 6439668
#> Projected CRS: RGF93 v1 / Lambert-93

Explore and display

  1. Plot municpalities with sf. The filling color must be blue and borders must be thick.
plot(st_geometry(mun), col = "lightblue", lwd = 2)

  1. Plot municipalities with mapsf. The filling color must be lightgreen and borders must be white and thin.
# install.packages("mapsf")
library(mapsf)
mf_map(mun, col = "lightgreen", border = "#ffffff", lwd = .5)

  1. Add the restaurants to the municpalities map.
rest <- st_read("data/lot.gpkg", layer = "restaurants")
#> Reading layer `restaurants' from data source 
#>   `/home/tim/Documents/prj/cartography_with_r/data/lot.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 694 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 542250 ymin: 6350950 xmax: 634198.4 ymax: 6437250
#> Projected CRS: RGF93 v1 / Lambert-93
mf_map(mun, col = "lightgreen", border = "#ffffff", lwd = .5)
mf_map(rest, pch = 21, col = "red", cex = 1, add = TRUE)

  1. How many municipalities are there?
# dim(mun)[1]
nrow(mun)
#> [1] 313
# length(unique(mun$INSEE_COM))
  1. Which commune has the largest population?
# mun[order(mun$POPULATION, decreasing = T), 2][1, ]
mun[mun$POPULATION == max(mun$POPULATION), 2]
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 570470.4 ymin: 6368056 xmax: 581068.4 ymax: 6380468
#> Projected CRS: RGF93 v1 / Lambert-93
#>    NOM_COM                           geom
#> 39  Cahors MULTIPOLYGON (((580994.5 63...

Attribute selection and join

  1. Import the com.csv file.
mun_df <- read.csv("data/com.csv")
  1. Join this dataset to the municipality layer.
mun_final <- merge(
  x = mun, # the sf object
  y = mun_df, # the data.frame
  by = "INSEE_COM", # identifier in x and y
  all.x = TRUE # keep all lines
)
  1. Extract municipalities with more than 500 active workers and whose share of active workers in the total population is greater than 30%.
mun_sel <- mun_final[mun_final$ACT > 500 & mun_final$SACT > 30, ]
  1. Plot a map with all municipalities in gray and selected municipalities in red.
mf_map(mun_final, col = "grey40", border = 'lightblue')
mf_map(mun_sel, col = "red", add = TRUE)

Geoprocessing

Spatial selection and join

  1. Import the layers of municipalities and restaurants.
rest <- st_read("data/lot.gpkg", layer = "restaurants")
#> Reading layer `restaurants' from data source 
#>   `/home/tim/Documents/prj/cartography_with_r/data/lot.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 694 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 542250 ymin: 6350950 xmax: 634198.4 ymax: 6437250
#> Projected CRS: RGF93 v1 / Lambert-93
mun  <- st_read("data/lot.gpkg", layer = "communes")
#> Reading layer `communes' from data source 
#>   `/home/tim/Documents/prj/cartography_with_r/data/lot.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 313 features and 12 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 539668.5 ymin: 6346290 xmax: 637380.9 ymax: 6439668
#> Projected CRS: RGF93 v1 / Lambert-93
  1. Perform a spatial join to find the name and identifier of the municipality where each restaurant is located.
resto <- st_join(x = rest, y = mun[, "INSEE_COM"], 
                 join = st_intersects,
                 left = FALSE)
resto
#> Simple feature collection with 694 features and 3 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 542250 ymin: 6350950 xmax: 634198.4 ymax: 6437250
#> Projected CRS: RGF93 v1 / Lambert-93
#> First 10 features:
#>      id                       nom INSEE_COM                     geom
#> 1   212              forest_cobra     46154 POINT (611092.9 6365662)
#> 2  2327  meteoropathologic_marlin     46003   POINT (596050 6414450)
#> 3  2328            pumice_unicorn     46003   POINT (596350 6414950)
#> 4  2329  indifferent_nakedmolerat     46003   POINT (595750 6417150)
#> 5  2330       pedagogic_eskimodog     46004   POINT (613350 6405150)
#> 6  2331  painstaking_alaskanhusky     46005 POINT (556763.5 6377394)
#> 7  2332 taxonomical_archaeopteryx     46006   POINT (573350 6411050)
#> 8  2333       spotless_hornedtoad     46008   POINT (560650 6389150)
#> 9  2334            dramatic_hyrax     46009   POINT (610450 6398050)
#> 10 2335     carefree_mountainlion     46009 POINT (610640.9 6397839)

Geometrical operations

  1. Compute the number of restaurants per municipality.
mun$n <- lengths(st_intersects(mun, resto, sparse = T))
########### OR
# resto <- st_join(x = rest, y = mun[, "INSEE_COM"], 
#                  join = st_intersects,
#                  left = FALSE)
# rest_by_mun <- aggregate(x = list(n = resto$id), 
#                          by = list(INSEE_COM = resto$INSEE_COM), 
#                          FUN = length)
# mun <- merge(mun, rest_by_mun, "INSEE_COM", all.x = TRUE)
# mun[is.na(mun$n), "n"] <- 0
  1. Which municipalities have more than 10 restaurants and fewer than 1000 inhabitants?
mun_sel <- mun[mun$n > 10 & mun$POPULATION < 1000, ]
  1. Create a map showing all the municipalities in grey and the municipalities selected above in red.
mf_map(mun_final, col = "grey40", border = 'lightblue')
mf_map(mun_sel, col = "red", add = TRUE)

Cartography

The Cartographic Language

How to represent the following variables:

  • A municipal population
    Proportional Symbols
  • A median age by region
    Choropleth map, ordered colors
  • A population growth rate
    Choropleth map, ordered colors
  • The administrative status of municipalities
    Typologie map, unordered colors
  • Life expectancy per country
    Choropleth map, ordered colors

mapsf

Map Types

  1. Import the com.csv file.
mun_df <- read.csv("data/com.csv")
  1. Join this dataset to the municipality layer.
mun  <- st_read("data/lot.gpkg", layer = "communes")
#> Reading layer `communes' from data source 
#>   `/home/tim/Documents/prj/cartography_with_r/data/lot.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 313 features and 12 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 539668.5 ymin: 6346290 xmax: 637380.9 ymax: 6439668
#> Projected CRS: RGF93 v1 / Lambert-93
mun <- merge(
  x = mun, # the sf object
  y = mun_df, # the data.frame
  by = "INSEE_COM", # identifier in x and y
  all.x = TRUE # keep all lines
)
  1. Create a map of the working population.
library(mapsf)
mf_map(mun, col = "ivory4", border = "ivory" , lwd = .3)
mf_map(
  x = mun  ,
  var = "ACT",
  type = "prop",
  inches = .35,
  symbol = "square",
  border = "white",
  lwd = .75,
  col = "#881094",
  leg_title = "Working Population",
  leg_pos = c(622700, 6372800)
)
mf_title("Working Population in the Lot Region")

  1. Create a map of the share of the working population in the total population.
mf_distr(mun$SACT)

summary(mun$SACT)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   7.812  23.104  26.754  27.265  31.537  52.083
# Classification method based mean and standard deviation
bks <- mf_get_breaks(mun$SACT, breaks = "msd", central = FALSE)
# Double color palette
cols <- mf_get_pal(c(3,4), palette = c("Mint", "Burg"))

mf_map(
  x = mun,
  var = "SACT",
  type = "choro",
  breaks = bks,
  pal = cols,
  border = "ivory4",
  lwd = .2,
  leg_title = "%",
  leg_val_rnd = 0
)
mf_title("Share of the working population in the total poulation")

Map Layout

  1. Create a map showing the working population in the industrial sector.
  2. Add the necessary layout elements.
  3. Make the map more intelligible, more explicit.
  4. Export the map in PNG format, 1000 pixels wide.

Workers + Workers in the industrial sector

dep <- st_read("data/lot.gpkg", layer = "departements", quiet = TRUE)
# Workers  + Workers in the industrial sector
par(mfrow = c(1, 2))

mf_map(x = mun, border = "white", lwd = .2, add = F)
mf_map(x = dep, lwd = 1, col = NA, add = TRUE, lend = 0)
mf_map(x = mun, var = "ACT", 
       type = "prop",
       inches=.2,
       leg_frame = T,
       leg_title = "N.")
mf_title("Working Population", tab = FALSE)
mf_credits(
  paste0("Admin Express COG Carto 3.0, IGN - 2021 & ",
         "BD CARTO® 4.0, IGN - 2021 ; Recensements harmonisés - ",
         "Séries départementales et munmunales, INSEE - 2020\n",
         "Auteurs : T. Giraud, 2025"), bg ="#ffffff80")
mf_map(x = mun, border = "white", lwd = .2, add = F)
mf_map(x = dep, lwd = 1, col = NA, add = TRUE, lend = 0)
mf_map(x = mun, var = "IND", 
       type = "prop",
       inches=.2,
       leg_frame = T,
       leg_title = "N.")
#> 50 '0' values are not plotted on the map.
mf_title("Working Population in the Industrial Sector", tab = FALSE)
mf_scale(5)
mf_arrow(pos = "topright")

Share of workers in the industrial secto

mf_distr(mun$SACT_IND)

summary(mun$SACT_IND)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   0.000   7.407  14.286  16.282  23.077 100.000

2 municipalities have a 100% share of workers in the industrial sector. Both have less than 15 workers. They are outliers.

# Remove outliers
mun_sel <- mun[mun$ACT > 15, ]
mf_distr(mun_sel$SACT_IND)

summary(mun_sel$SACT_IND)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   0.000   7.408  14.286  15.783  22.466  60.000
bks <- mf_get_breaks(mun_sel$SACT_IND, nbreaks = 5, breaks = "quantile")
mf_map(x = mun,
       var = "SACT_IND",
       type = "choro",
       breaks = bks,        
       leg_val_rnd = 0,      
       pal = "Red-Yellow",  
       leg_title = "%",
       add = FALSE,
       col_na = "grey",
       leg_no_data = "less than 15 workers")
#> 2 values are outside the class limits and will be classified as 'NA'.
mf_title("Share of workers in the industrial sector")
mf_scale(5)
mf_arrow(pos = "topright")
mf_credits(paste0("Auteurs : T. Giraud, 2025\n","Admin Express COG Carto 3.0, IGN - 2021 & ",
                  "BD CARTO® 4.0, IGN - 2021 ;\nRecensements harmonisés - ",
                  "Séries départementales et munmunales, INSEE - 2020"))

Final Map

mf_export(x = mun, filename = "img/map.png", width = 1000, res= 120)
mf_map(x = mun, border = "white", lwd = .2, add = T)
mf_map(x = dep, lwd = 1, col = NA, add = TRUE, lend = 0)
mf_map(mun, c("ACT", "SACT_IND"), "prop_choro",
       breaks = bks,
       pal = "Red-Yellow",
       inches = .4,
       border = "white", lwd = .7,
       leg_val_rnd =  c(0,1),
       leg_pos = c(538000,6442000, 538000, 6424000), 
       leg_title = c("N. workers",
                     "Workers in the\nIndustrial Sector\n(%)"),
       col_na = "grey",
       leg_no_data = "Less than 15 workers")
#> 2 values are outside the class limits and will be classified as 'NA'.
# Ajout d'annotations
mf_annotation(x = mun[mun$NOM_COM=="Biars-sur-Cère",], txt = "Andros\n(jam factory)",
             halo = TRUE, cex = 1)
mf_annotation(x = mun[mun$NOM_COM=="Figeac",], txt = "Aerospace\nIndustry",
             pos = "bottomright", halo = TRUE, cex = 1)
mf_annotation(x = mun[mun$NOM_COM=="Gramat",], txt = "La Quercynoise\n(duck meat factory)",
              pos = "topleft", halo = TRUE, cex = 1)

mf_title("MECANIC VALLEE")

# inset->
mf_inset_on(fig = c(.8,0.98,0.1,0.3))
mf_map(dep, lwd = .1)
mf_map(mun, border = NA, add = TRUE, col = "#D7FF68")
box(lwd = .5)
mf_inset_off()
# <- inset
mf_scale(5)
mf_arrow("topright")
mf_credits(paste0("Auteurs : T. Giraud, 2025\n","Admin Express COG Carto 3.0, IGN - 2021 & ",
                  "BD CARTO® 4.0, IGN - 2021 ;\nRecensements harmonisés - ",
                  "Séries départementales et munmunales, INSEE - 2020"))
dev.off()
#> png 
#>   2